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A fast high-order method to calculate wakefield forces in an electron beam 

Ji Qiang, Chad Mitchell, Robert D. Ryne 
Lawrence Berkeley National Laboratory, 
1 Cyclotron Road, Berkeley, CA 94720 

In this paper we report on a high-order fast method to numerically calculate wakefield 



Vh ' forces in an electron beam given a wake function model. This method is based on a Newton- 

Cotes quadrature rule for integral approximation and an FFT method for discrete summation 



that results in an 0{Nlog{N)) computational cost, where N is the number of grid points. 
^^ , Using the Simpson quadrature rule with an accuracy of 0{h'^), where h is the grid size, we 

^ ■ present numerical calculation of the wakefields from a resonator wake function model and 

C*H' from a one-dimensional coherent synchrotron radiation (CSR) wake model. Besides the fast 
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o 



speed and high numerical accuracy, the calculation using the direct line density instead of the 
first derivative of the line density avoids numerical filtering of the electron density function 



C/3 

O ■ for computing the CSR wakefield force. 

Ph; I. INTRODUCTION 

^ , The wakefield forces from the interaction of electrons with external materials and from the 

in 

(^ ■ interactions of electrons inside the beam through radiation have significant effects on the electron 

m : 

1/^ , beam quality for next generation light sources. Those wakefields can cause the loss of beam energy, 



increase the beam energy spread, increase beam emittance through the magnetic bunch compressor, 
and drive microbunching instability through the linac. Including the wakefield effects in a beam 
dynamics tracking code is crucial for designing a beam delivery system for next generation light 
p\ . sources and has been included in a number of tracking codes [11-14]. 

The numerical calculation of the wakefield forces on electrons involves calculating a convolution 
between the wake function and the charge density of the beam. The direct numerical calculation 
of the convolution for the wakefield has a computational cost scaling as 0{N'^), where N is the 
number of grid points in the computational domain. Fortunately, such a discretized convolution 
summation on a uniform grid can be calculated using a cyclic summation on a doubled computa- 
tional domain using an FFT based method 5|-l8(. This reduces the computational cost from 0{N'^) 
to 0{Nlog{N)). In our previous study, we have developed a high-order FFT based method to nu- 
merically calculate the convolution with a smooth kernel [9|. In this paper, we extend the previous 



study by applying a high-order Simpson quadrature rule to the wakefield convolution with finite 
discontinuity of the kernel. We then apply this method to the calculation of the wakefields from 
an external resonator wake function (an approximation to the AC resistance wall wake function) 
and from a steady-state CSR wake function. 

The organization of the paper is as follows: after a brief introduction, we present the com- 
putational method in Section 2, application examples to calculate the wakefield forces from an 
external resonator wake function and from a steady-state CSR wake function in Section 3, and 
final summary in Section 4. 

II. COMPUTATIONAL METHOD 

We consider the longitudinal wakefield on an electron at longitudinal location z inside an electron 
beam can be written as 

/oo 
W{z - z')p{z)dz' (1) 

where W is the wake function from the interaction of electrons with external materials such as 
resistive wall wake or the interaction among electrons themselves such as the coherent synchrotron 
radiation (CSR) inside a bending magnet, and p is the beam line density distribution. For the 
external wake function, the variable z stands for ct, where t is the time of flight difference with 
respect to the head of the beam, and c is the speed of light in vacuum. For the CSR wake function, 
the variable z stands for the longitudinal spatial position. The above convolution integral can be 
written as the summation of M equal subinterval integrals: 

^(^) = E/ '^'W{z-z')p{z')dz' (2) 

where z'^ = {i — l){z/M),i = 1, • • • ,M. For the integral between z'^ and z'^^^, the closed Newton- 
Cotes formula of degree k at k + 1 equally spaced points can be written as 

r^' W{z - z')p{z')dz' « Y. ^^-^(^ - ^'M^'j) (3) 

where z', = z'^+jdz, dz = z/{kM), and the weight Wj is associated with integral of the jth Lagrange 
basis polynomial. For k = 1, this is known as the trapezoidal quadrature rule; for k = 2, this is the 



Simpson rule; for A; = 3, Simpson's 3/8 rule; for k = 4, Boole's rule, etc [lOl,!!!!- Using the extended 
three-point Simpson rule on a uniform grid, the above convolution integral can be approximated 
as 

Eiz) = E{z) + 0{h^) (4) 



where h is the grid size, and E{z) is the numerical approximation of the convolution integral. For 
a grid point Zj, if i is an odd number, the Simpson quadrature rule yields 

(i-l)/2 (i-l)/2 



1 
3' 



E(z,.) = ^M/i+4 E /2,+2 J2 hj-i + fi. 
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if i is an even number, it yields 



(i-2)/2 i/2 

^(^0 = ^/i(3/i/2 + 5/2/2 + 2 5] /2. + 4^/2,_i + /,: 
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where the kernel function / is given by 



(5) 



(6) 



f{z,z') = W{z-z')p{z') (7) 

and fj{zi) = f{zi,z'j) with z' = (j — l)/i. The above numerical summation can be rewritten as 



1 ^ 1 

i?(z,) = -hJ2w{z,-z,)p{^,) + -hW{0)p{z,) 



with 
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for odd grid number i, and 
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(8) 



(9) 



(10) 



(11) 



The direct calculation of the above summation requires 0{N) operations for a single point Zi, 
where z = 1, • • • , iV. To obtain the convolution for all N points on the grid, the total computational 
cost will be 0{N'^). Fortunately, the above convolution summation can be replaced by a cyclic 
summation in the double-gridded computational domain: 

2N 



1 



Ecizi) = i:h^Wc{zi- Zj)pc 



(12) 



i=i 



where i = 1, • • • , 2iV, j = 1, • • • , 2iV and 

Pc{zj) = { (13) 

: N <j <2N 

W(zk) : 1< A; < iV; 
Ty(-Z2iv+2-fc) : N<k<2N; 
where the new density function pdz) and the new wake function Wdz) have the properties: 

Pciz) = pciz + 2iL + h)) (15) 

Wc{z) = Wdz + 2{L + h)) (16) 

where L = (N — l)h. From above definition, one can show that the cychc summation gives the 
same value as the convolution summation within the original domain, i.e. 

E{zi) = E^Zi) for 1 = 1, N. (17) 

In the cyclic summation, the kernel is a discrete periodic function in the doubled computational 
domain. This cyclic summation can be calculated using the FFT method, i.e. 

E^Zi) = InvFFT{Ed^)) (18) 

Here, InvFFT denotes the inverse FFT of the function Ec{uj) that is given by 

E^oo) = Wciu) pdoo) (19) 

where Wc{uj) and Pc{oj) denote the forward FFT of the function Wc and the pc respectively. The 
computational operations required to calculate cyclic summation using the above FFT method is 
oiO{Nlog(N). 

III. APPLICATION EXAMPLES 

As an illustration, we use the above method to calculate the wakefields inside an electron beam 
from a resonator wake function and from the CSR wake function. The resonator wake function is 
given as [Ij] 

W{s) = 4exp(-(s/so)/r)cos((8/r)i/4(s/so)); for s>0 (20) 

where s = z — z' is the separation between the observation point and the source point, sq = 
{2a/[ZQa))^''^ is the characteristic length of the wake, a is the radius of the pipe, a is the conduc- 
tivity of the pipe material, Zq is the vacuum impedance, F = ct/sq is the normalized relaxation 




FIG. 1: The wakefield force as a function of bunch length coordinate in an electron beam inside a Copper 
pipe with resonator wake function. 
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FIG. 2: The errors at the center of the electron beam as a function grid points. The red line is from the 
Simpson quadrature rule and the green line is from the trapzoidal rule. 



time, and r is the relaxation time. This wake function provides a good approximation to the AC 
resistive wall wake function at large T J12l |. 

Figure [1] shows the normalized wakefield in an electron beam moving inside a Copper pipe 
(r = 1) computed using the above algorithm. The beam has an rms bunch length of I.Ssq and a 
Gaussian current distribution. The electron beam loses energy to the wall due to the resistance 
of the conducting pipe. To verify the accuracy of the above algorithm, Figure [2] shows the errors 
of the calculated wakefield at the center of the bunch as a function of the number of grid points. 
As a comparison, we also give the errors using the trapezoidal quadrature rule. As expected, the 
above Simpson rule method converges much faster than the trapezoidal rule based method. The 
accuracy of the Simpson rule calculated wakefield increases by more than an order of magnitude 
as the number of grid points doubles. 

Another application example of above numerical method is to calculate the wakefield force from 
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FIG. 3: The GSR wakefield as a function of bunch length coordinate in an electron beam. 

the one-dimensional steady state CSR wake. The one-dimensional CSR wake function is given 
as []J 

1/6 - uVl8 - u'^/64 



W{U[S)) - 4^g ^2 i 2(1+ ^2/4)3 + (1+^2/4)3(1 + ^2/12)2 



-} 



and 



(s) = (\/64 + 144s2 + l2s)^/3 - (V64 + 144^2 - 12s)^ 



/3 



(21) 



(22) 



where s = ^ is the normalized longitudinal separation between the source point and the obser- 
vation point, R is the radius of the bending magnet, and 7 is the relativistic factor of the electron 
beam. The direct application of the above method with different number of grid points is shown 
in Fig. [3] for a 1 nC Gaussian beam with 10 micron rms bunch length at 100 MeV. It is seen that 
a few thousand grid points are needed in order to obtain a numerically converged wakefield. 
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FIG. 4: The steady-state GSR wake function as a function of separation between the observation point and 
the source point. 



This large number of grid points is needed due to the sharp variation of the CSR wake function 
around s = as shown in Fig. HI It is seen that outside this short range (about a few micron in this 



example) around zero, the CSR wake function is much smoother and can be approximated by an 
asymptotic function s~^'^. Recently, an integrated Green function method with linear basis was 
proposed to treat the CSR wake function including the short range effect IJ]. Here, we propose 
a modified Simpson rule method by separating the range of the original CSR wake function into 
two regions: one small region around the zero covering the short range CSR interaction, and the 
remainning region covering the long range interaction. The integral for the CSR wakefield can be 

written as: 

rz—5z 



E{z) 



W{z-z')p(z')dz' + I W{z-z')p{z)dz 

-5z 



(23) 



where 5z is a small distance around the zero to separate the long-range CSR interaction from the 
short-range interaction. The first integral in the above equation can be calculated using a Simpson 
quadrature rule as shown in the previous section. Using the parameters in the above example with 
5z = 3.2 um and the Simpson quadrature rule, we calculated the errors at the center of the bunch 
as a function of the number of grid points for the first integral in Eq. [23l The results are shown 
in Fig. [5l As expected, the error of the integral decreases by more than an order of magnitude as 
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FIG. 5: The errors at the center of the bunch as a function of grid points for the first integral in Eq. 
using the Simpson rule. 



the number of the grid points doubles. 

The second integral in the Eq. [23] will be calculated using an integrated Green's function method 
with constant current density basis functions given the fact that the CSR wake function varies 
sharply while the line density changes slowly inside such a short-range interval. At each grid 
point Zi, the interval [zi — 6z,Zi] is subdivided into Nmt number of sub-intervals with a length 
hint = Sz/(Nint — !)• This integral can now be rewritten as: 

W{zi - z')p{z')dz' = / W{zi - z')p{z')dz' + J2 ^(^» - z')p{z')dz' 

Sz Jzi ■2J Zj—hint/'i 



+ / W{zi - z')p{z')dz' 



(24) 



where zi = Zi — bz. Since the bz is small (0(/U7ti-)), only a small number of sub-intervals are needed. 
Assuming that the current density is constant within each sub-interval, the above integral can be 
approximated as: 

Zl+hint/'i , , ^Z^ pZj+hint/2 



W{zi - z')p{z')dz' = p{zi) 

Zi—Sz J Z\ 



W{zi - z')dz' + ^ p{zj) I W{zi - z)dz 



+p{zi) / W{zi - z')dz' 



(25) 



The integral with respect to the wake function in above equation can be done numerically if an 
analytical solution is not available. Fortunately, in the case of the steady-state CSR wake function, 
there exists an analytical solution for the integral of the wake function [l5|, i.e. 



Zj l^int / ^ 



W{z-Z')dz' = [U{A(Z -Zj + hint/2)) -v{A{z-Zj-hi^^/2))]/ A 



where 
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^ :(f^l/3+J^-l/3)+ 2 



:(f^2/3 _ f^-2/3)| 



(26) 



(27) 



and Q = u + \/v? + l, and A = 3"f'^/{2R). Using this modified Simpson rule, we recalculated the 
CSR Wakefield inside the electron beam using 10 sub-intervals inside the short range around and 
127 grid points for the rest of the region. The results are shown in Fig. [6] together with the direct 
Simpson rule calculation using 8191 grid points. It is seen that the wakefield from the modified 
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FIG. 6: The CSR wakefield as a function of bunch length coordinate in an electron beam with the modified 
Simpson rule (137 grid points) and the direct Simpson rule (8191 grid points). 



Simpson rule agrees with the high resolution direct Simpson rule very well even with only 137 grid 
points. 



IV. SUMMARY 

In this paper, we present a fast high-order method for numerical calculation of the wakefield 
on electrons inside an electron beam. Using the Simpson quadrature rule, we present examples 
to calculate the wakefields from the resonator wake function and from the CSR wake function. 
These examples can have 0{h^) accuracy with a 0{Nlog(N)) computational cost. Besides the fast 
speed and potential high numerical accuracy, the above CSR wakefield calculation also uses the 
direct line density instead of the first derivative of the line density with numerical filtering of the 
density as suggested in Reference [16(]. This paper gives application examples to the calculation 
of the longitudinal wakefield for an electron beam. The same method can also be applied to the 
calculation of the transverse wakefields inside an electron beam using a modified density function. 
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